﻿using SoftwareConsulting.BI11.NumericalMethods;

namespace SoftwareConsulting.BI11.NumericalMethods.MatrixAlgebra
{
    /// <summary>
    /// Класс симметричной матрицы
    /// </summary>
    public class SymmetricMatrix : Matrix
    {
        private const int _lowerUpperPermutationCRLCriticalDimension = 36;

        /// <summary>
        /// Конструктор класса
        /// </summary>
        /// <param name="matrixComponents">элементы матрицы</param>
        internal protected SymmetricMatrix(double[,] matrixComponents) : 
            base(matrixComponents)
        {
        }

        /// <summary>
        /// Конструктор класса
        /// </summary>
        /// <param name="matrixDimensionality">размерность матрицы</param>
        public SymmetricMatrix(int matrixDimensionality) : 
            base(matrixDimensionality, matrixDimensionality)
        {
        }

        /// <summary>
        /// Конструктор класса
        /// </summary>
        /// <param name="rowsNumber">число строк</param>
        /// <param name="columnsNumber">число столбцов</param>
        public SymmetricMatrix(int rowsNumber, int columnsNumber) : 
            base(rowsNumber, columnsNumber)
        {
        }

        /// <summary>
        /// Возвращает симметричную матрицу, полученную сложением двух симметричных матриц
        /// </summary>
        /// <param name="B"></param>
        /// <param name="A"></param>
        /// <returns></returns>
        public static SymmetricMatrix operator+(SymmetricMatrix B, SymmetricMatrix A)
        {
            return new SymmetricMatrix(B.GetResultOfComponentsAddition(A));
        }

        /// <summary>
        /// Возвращает обратную матрицу, используя алгоритм CLR
        /// </summary>
        /// <returns></returns>
        private SymmetricMatrix GetInversedMatrixUsingCLRAlgorithm()
        {
            if (this.RowsCount == 1)
                return GetInversed1x1Matrix();
            if (this.RowsCount == 2)
                return GetInversed2x2Matrix();
            Matrix[] splitMatrices = GetSplittedMatrixes();
            SymmetricMatrix b1 = (SymmetricMatrix)splitMatrices[0].GetInversedMatrix();
            Matrix cb1 = splitMatrices[2].GetResultOfMultiplication(b1);
            SymmetricMatrix cb1cT = new SymmetricMatrix(cb1.GetResultOfComponentsTransposedMultiplication(splitMatrices[2]));
            splitMatrices[1] = ((SymmetricMatrix)splitMatrices[1]).GetResultOfSubtraction(cb1cT).GetInversedMatrix();
            splitMatrices[2] = splitMatrices[1].GetResultOfMultiplication(cb1);
            splitMatrices[0] = b1.GetResultOfAddition(new SymmetricMatrix(cb1.GetResultOfComponentsSelfTransposedMultiplication(splitMatrices[2])));
            return SymmetricMatrix.Join(splitMatrices);
        }

        /// <summary>
        /// Возвращает симметричную матрицу по указанным элементам
        /// </summary>
        /// <param name="matrixComponents"></param>
        /// <returns></returns>
        public static SymmetricMatrix FromComponents(double[,] matrixComponents)
        {
            if (matrixComponents.GetLength(0) != matrixComponents.GetLength(1))
                throw new BIMatrixNotSquareException(matrixComponents.GetLength(0), matrixComponents.GetLength(1), null);
            for (int i = 0; i < matrixComponents.GetLength(0); i++)
            {
                for (int j = 0; j < i; j++)
                {
                    if (matrixComponents[i, j] != matrixComponents[j, i])
                        throw new BIMatrixNonSymmetricComponentsException(i, j, matrixComponents[i,j], matrixComponents[j, i], null);
                }
            }
            return new SymmetricMatrix(matrixComponents);
        }

        /// <summary>
        /// Возвращает единичную симметричную матрицу указанной размерности
        /// </summary>
        /// <param name="matrixDimensionality"></param>
        /// <returns></returns>
        public static SymmetricMatrix IdentityMatrix(int matrixDimensionality)
        {
            double[,] a = new double[matrixDimensionality, matrixDimensionality];
            for (int i = 0; i < matrixDimensionality; i++)
            {
                a[i, i] = 1;
            }
            return new SymmetricMatrix(a);
        }

        /// <summary>
        /// Возвращает обратную матрицу
        /// </summary>
        /// <returns></returns>
        public override Matrix GetInversedMatrix()
        {
            return this.RowsCount < _lowerUpperPermutationCRLCriticalDimension
                        ? new SymmetricMatrix(new LowerUpperPermutationDecomposition(this).GetInversedMatrixComponents())
                        : GetInversedMatrixUsingCLRAlgorithm();
        }

        /// <summary>
        /// Возвращает обратную матрицу
        /// </summary>
        /// <returns></returns>
        /// <remarks>Применимо только для случая матрицы 1x1 !</remarks>
        private SymmetricMatrix GetInversed1x1Matrix()
        {
            double[,] newComponents = new double[1, 1];
            newComponents[0, 0] = 1 / Components[0, 0];
            return new SymmetricMatrix(newComponents);
        }

        /// <summary>
        /// Возвращает обратную матрицу
        /// </summary>
        /// <returns></returns>
        /// <remarks>Применимо только для случая матрицы 2x2 !</remarks>
        private SymmetricMatrix GetInversed2x2Matrix()
        {
            double[,] newComponents = new double[2, 2];
            double inverseDeterminant = 1 / (Components[0, 0] * Components[1, 1]
                                        - Components[0, 1] * Components[1, 0]);
            newComponents[0, 0] = inverseDeterminant * Components[1, 1];
            newComponents[1, 1] = inverseDeterminant * Components[0, 0];
            newComponents[0, 1] = newComponents[1, 0] = -inverseDeterminant * Components[1, 0];
            return new SymmetricMatrix(newComponents);
        }        

        /// <summary>
        /// Возвращает симметричную матрицу, полученную умножением симметричной матрицы на коэффициент
        /// </summary>
        /// <param name="A"></param>
        /// <param name="C"></param>
        /// <returns></returns>
        public static Matrix operator *(SymmetricMatrix A, double C)
        {
            return new SymmetricMatrix(A.GetResultOfComponentsMultiplication(C));
        }

        /// <summary>
        /// Возвращает симметричную матрицу, полученную перемножением двух симметричных матриц
        /// </summary>
        /// <param name="B"></param>
        /// <param name="A"></param>
        /// <returns></returns>
        public static SymmetricMatrix operator*(SymmetricMatrix B, SymmetricMatrix A)
        {
            return new SymmetricMatrix(B.GetResultOfComponentsMultiplication(A));
        }

        /// <summary>
        /// Возвращает симметричную матрицу, полученную умножением на симметричную матрицу, транспонированную относительно указанной
        /// </summary>
        /// <param name="secondOriginalMatrix"></param>
        /// <returns></returns>
        public SymmetricMatrix GetResultOfTransposedMultiplication(SymmetricMatrix secondOriginalMatrix)
        {
            if (secondOriginalMatrix.ColumnsCount != this.ColumnsCount)
                throw new BIMatrixCantMultiplyTransposedMatrixException(this.ColumnsCount, secondOriginalMatrix.ColumnsCount, null);
            return new SymmetricMatrix(GetResultOfComponentsTransposedMultiplication(secondOriginalMatrix));
        }

        /// <summary>
        /// Возвращает симметричную матрицу, полученную сложением с указанной симметричной матрицей
        /// </summary>
        /// <param name="secondMatrix"></param>
        /// <returns></returns>
        protected SymmetricMatrix GetResultOfAddition(SymmetricMatrix secondMatrix)
        {
            return new SymmetricMatrix(GetResultOfComponentsAddition(secondMatrix));
        }	

        /// <summary>
        /// Возвращает симметричную матрицу, полученную умножением на указанную симметричную матрицу
        /// </summary>
        /// <param name="secondMatrix"></param>
        /// <returns></returns>
        protected SymmetricMatrix GetResultOfMultiplication(SymmetricMatrix secondMatrix)
        {
            return new SymmetricMatrix(GetResultOfComponentsMultiplication(secondMatrix));
        }

        /// <summary>
        /// Возвращает симметричную матрицу, полученную вычитанием указанной симметричной матрицы
        /// </summary>
        /// <param name="secondMatrix"></param>
        /// <returns></returns>
        protected SymmetricMatrix GetResultOfSubtraction(SymmetricMatrix secondMatrix)
        {
            return new SymmetricMatrix(GetResultOfComponentsSubtraction(secondMatrix));
        }

        /// <summary>
        /// Возвращает разбивку на 3 матрицы примерно одинакового размера
        /// </summary>
        /// <returns></returns>
        private Matrix[] GetSplittedMatrixes()
        {
            int n = this.RowsCount;
            int p = MathExtended.LargestPowerOf2SmallerThan(n);
            int q = n - p;
            double[,] a = new double[p, p];
            double[,] b = new double[q, q];
            double[,] c = new double[q, p];
            for (int i = 0; i < p; i++)
            {
                for (int j = 0; j < p; j++)
                    a[i, j] = Components[i, j];
                for (int j = p; j < n; j++)
                    c[j - p, i] = Components[i, j];
            }
            for (int i = p; i < n; i++)
            {
                for (int j = p; j < n; j++)
                    b[i - p, j - p] = Components[i, j];
            }
            return new Matrix[]
                    { new SymmetricMatrix(a), new SymmetricMatrix(b), new Matrix(c) };
        }

        /// <summary>
        /// Строит матрицу из трех частей (обратно к разбирению)
        /// </summary>
        /// <param name="parts"></param>
        /// <returns></returns>
        private static SymmetricMatrix Join(Matrix[] parts)
        {
            int p = parts[0].RowsCount;
            int n = p + parts[1].RowsCount;
            double[,] newComponents = new double[n, n];
            for (int i = 0; i < p; i++)
            {
                for (int j = 0; j < p; j++)
                    newComponents[i, j] = parts[0][i, j];
                for (int j = p; j < n; j++)
                    newComponents[i, j] = newComponents[j, i] = -parts[2][j - p, i];
            }
            for (int i = p; i < n; i++)
            {
                for (int j = p; j < n; j++)
                    newComponents[i, j] = parts[1].Components[i - p, j - p];
            }
            return new SymmetricMatrix(newComponents);
        }        

        /// <summary>
        /// Возвращает разность двух симметричных матриц
        /// </summary>
        /// <param name="B"></param>
        /// <param name="A"></param>
        /// <returns></returns>
        public static SymmetricMatrix operator-(SymmetricMatrix B, SymmetricMatrix A)
        {
            return new SymmetricMatrix(B.GetResultOfComponentsSubtraction(A));
        }

        /// <summary>
        /// Возвращает транспонированную симметричную матрицу
        /// </summary>
        /// <returns></returns>
        /// <remarks>т.е. саму себя</remarks>
        public override Matrix GetTransposedMatrix()
        {
            return this;
        }

        /// <summary>
        /// Возвращает симметричную матрицу, полученную умножением своей транспонированной матрицы на указанную матрицу
        /// </summary>
        /// <param name="secondMatrix"></param>
        /// <returns></returns>
        public SymmetricMatrix GetResultOfSelfTransposedMultiplication(SymmetricMatrix secondMatrix)
        {
            if (secondMatrix.RowsCount != this.RowsCount)
                throw new BIMatrixCantMultiplySelfTransposedMatrixException(this.RowsCount, secondMatrix.RowsCount, null);                        
            return new SymmetricMatrix(GetResultOfComponentsSelfTransposedMultiplication(secondMatrix));
        }
    }
}
